rm(list = ls())

library(tidyverse)
library(stargazer)
library(lfe)
library(broom) #for tidy
library(lmtest) # For coeftest()
library(sandwich) # For sandwich()
library(estimatr)
library(texreg)
library(psych)
library(default)
library(sjmisc) #for row_sums
library(quantreg)
library(fixest)

options(tibble.print_min = 10)
source("functions.R")

load("./data_output/persons_years_sample.Rdata")

full_sample %>% summarise(mean(inc_labor == 0))
full_sample_pos = full_sample %>% filter(inc_labor != 0)

full_sample_pos = full_sample_pos %>%  mutate_at(vars(extro:cog), myScale)

###########################################################################
# Estimations -------------------------------------------------------------
###########################################################################


# Earnings log
trendsPlot(full_sample_pos, "inc_log ~ cog + noncog")
ggsave("/figures/panel_inclog=cog+noncog.pdf", height = 4, width = 7)

trendsPlot(full_sample_pos, "inc_log ~ extro + consci") + labs(title = "Panel B: Noncognitive Skills Without Cognitive Skills")
ggsave("/figures/panel_inclog=extro+consci.pdf", height = 4, width = 7)

trendsPlot(full_sample_pos, "inc_log ~ cog + extro + consci") + labs(title = "Panel A: Noncognitive Skills Decomposed")
ggsave("/figures/panel_inclog=cog+extro+consci.pdf", height = 4, width = 7)

# Extensive margin
trendsPlot(full_sample, "extensive ~ cog + noncog") +
  theme(legend.position = "none")
ggsave("/figures/panel_extensive=cog+noncog.pdf", height = 4, width = 5.5)

trendsPlot(full_sample, "extensive ~ extro + consci") +
  theme(legend.position = "none")
ggsave("/figures/panel_extensive=extro+consci.pdf", height = 4, width = 5.5)

trendsPlot(full_sample, "extensive ~ cog + extro + consci") +
  theme(legend.position = "none")
ggsave("/figures/panel_extensive=cog+extro+consci.pdf", height = 4, width = 5.5)
# Earnings log intensive margin
trendsPlot(full_sample %>% filter(extensive == T), "inc_log ~ cog + noncog")
ggsave("/figures/panel_inclog_intensive=cog+noncog.pdf", height = 4, width = 6.7)

trendsPlot(full_sample %>% filter(extensive == T), "inc_log ~ extro + consci")
ggsave("/figures/panel_inclog_intensive=extro+consci.pdf", height = 4, width = 7)

trendsPlot(full_sample %>% filter(extensive == T), "inc_log ~ cog + extro + consci")
ggsave("/figures/panel_inclog_intensive=cog+extro+consci.pdf", height = 4, width = 7)




# Earnings rank
trendsPlot(full_sample, "inc_rank ~ cog + noncog") +
  theme(legend.position = "none")
ggsave("/figures/panel_incrank=cog+noncog.pdf", height = 4, width = 5.5)

trendsPlot(full_sample, "inc_rank ~ extro + consci") +
  theme(legend.position = "none")
ggsave("/figures/panel_incrank=extro+consci.pdf", height = 4, width = 5.5)

trendsPlot(full_sample, "inc_rank ~ cog + extro + consci") +
  theme(legend.position = "none")
ggsave("/figures/panel_incrank=cog+extro+consci.pdf", height = 4, width = 5.5)
# Earnings proportional
trendsPlot(full_sample, "inc_prop ~ cog + noncog")
ggsave("/figures/panel_incprop=cog+noncog.pdf", height = 4, width = 6.7)

trendsPlot(full_sample, "inc_prop ~ extro + consci")
ggsave("/figures/panel_incprop=extro+consci.pdf", height = 4, width = 7)

trendsPlot(full_sample, "inc_prop ~ cog + extro + consci")
ggsave("/figures/panel_incprop=cog+extro+consci.pdf", height = 4, width = 7)



# Education
trendsPlot(full_sample_pos, "inc_log ~ cog + consci + extro + educ_y") + labs(title = "Panel A: Skill Premia Conditional on Education")
ggsave("/figures/panel_inclog=cog+extro+consci+educ_y.pdf", height = 4, width = 7)

trendsPlot(full_sample, "educ_y ~ cog + consci + extro") + labs(title = "Panel B: Selection into Education")
ggsave("/figures/panel_educ_y=cog+consci+extro.pdf", height = 4, width = 7)


# Robustness

trendsPlot(full_sample_pos %>% mutate(extro = extro_deming), "inc_log ~ cog + extro + consci") + 
  labs(title = "Panel A: Extraversion = Sociability + Activity + Leadership + Confidence")
ggsave("/figures/robust_panel_inclog=cog+extro_deming+consci.pdf", height = 4, width = 7)

trendsPlot(full_sample_pos %>% mutate(extro = extro_factor, consci = consci_factor), 
           "inc_log ~ cog + extro + consci") +
  labs(title = "Panel B: Conscientiousness = Dutifulness + Deliberation\nExtraversion = Sociability + Activity + Leadership + Confidence + Achievement")
ggsave("/figures/robust_panel_inclog=cog+extro_factor+consci_factor.pdf", height = 4, width = 7)

trendsPlot(full_sample_pos, "inc_log ~ cog + extro + consci + leader_conf")+
  labs(title = "Panel C: Added Composite =  Leadership + Confidence")
ggsave("/figures/robust_panel_inclog=cog+extro+consci+leader_conf.pdf", height = 4, width = 7.1)

trendsPlot(full_sample_pos, "inc_log ~ cog + extro + consci + validity + lie + fix")
ggsave("/figures/robust_panel_inclog=cog+extro+consci+validity+lie+fix.pdf", height = 4, width = 7)


# Interactions
trendsPlot(full_sample_pos, "inc_log ~ cog*noncog")
ggsave("/figures/panel_inclog=cogXnoncog.pdf", height = 4, width = 6.7)

trendsPlot(full_sample_pos, "inc_log ~ cog*extro*consci - cog:extro:consci") #dont estimate triple interaction
ggsave("/figures/panel_inclog=cogXextroXconsci.pdf", height = 4, width = 7)

trendsPlot(full_sample_pos %>% mutate(educ_y = scale(educ_y)), "inc_log ~ cog*educ_y + noncog*educ_y")
ggsave("/figures/panel_inclog=cogXnoncogXeduc_y.pdf", height = 4, width = 7)

trendsPlot(full_sample_pos %>% mutate(educ_y = scale(educ_y)), "inc_log ~ cog*educ_y + extro*educ_y + consci*educ_y")
ggsave("/figures/panel_inclog=cogXextroXconsciXeduc_y.pdf", height = 4, width = 7)



# Additional --------------------------------------------------------------

#trendsPlot(full_sample %>% filter(extensive == T), "inc_log ~ cog + consci + extro + educ_y")
#ggsave("/figures/panel_inclog_intensive=cog+extro+consci+educ_y.pdf", height = 4, width = 7)


# # Employment
# trendsPlot(full_sample, "emp_11 ~ cog + noncog")
# ggsave("/figures/panel_emp=cog+noncog.pdf", height = 4, width = 7)
# 
# trendsPlot(full_sample, "emp_11 ~ extro + consci")
# ggsave("/figures/panel_emp=extro+consci.pdf", height = 4, width = 7)
# 
# trendsPlot(full_sample, "emp_11 ~ cog + extro + consci")
# ggsave("/figures/panel_emp=cog+extro+consci.pdf", height = 4, width = 7)
# 
# trendsPlot(full_sample, "emp_11 ~ cog + consci + extro + educ_y")
# ggsave("/figures/panel_emp=cog+extro+consci+educ_y.pdf", height = 4, width = 7)

# trendsPlot(full_sample, "educ_y ~ cog + noncog")
# ggsave("/figures/panel_educ_y=cog+noncog.pdf", height = 4, width = 7)

# # Robustness
# trendsPlot(full_sample_pos %>% mutate(extro = extro_deming, consci = consci_factor), 
#            "inc_log ~ cog + extro + consci")
# ggsave("/figures/robust_panel_inclog=cog+extro_deming+consci_factor.pdf", height = 4, width = 7)
